1 



Constructing a Stochastic Model of Bumblebee Flights from 
Experimental Data 

Friedrich Lcnz^'*, Alcksci V. Chcchkin^-'^, Raincr Klagcs-'^ 

1 School of Mathematical Sciences, Queen Mary University of London, Mile End Road, 
London El 4NS, UK 

2 Max-Planck Institute for Physics of Complex Systems, Noethnitzer Str. 38, 01187 
Dresden, Germany 

3 Institute for Theoretical Physics, NSC KIPT, ul. Akademicheskaya 1,UA-61108 
Kharkov, Ukraine 

* E-mail: f.lenz@qmul.ac.uk 

Abstract 

The movement of organisms is subject to a multitude of influences of widely varying character: from 
the bio-mechanics of the individual, over the interaction with the complex environment many animals 
live in, to evolutionary pressure and energy constraints. As the number of factors is large, it is very 
hard to build comprehensive movement models. Even when movement patterns in simple environments 
are analysed, the organisms can display very complex behaviours. While for largely undirected motion 
or long observation times the dynamics can sometimes be described by isotropic random walks, usually 
the directional persistence due to a preference to move forward has to be accounted for, e.g., by a 
correlated random walk. In this paper we generalise these descriptions to a model in terms of stochastic 
differential equations of Langevin type, which we use to analyse experimental search flight data of foraging 
bumblebees. Using parameter estimates we discuss the differences and similarities to correlated random 
walks. From simulations we generate artificial bumblebee trajectories which we use as a validation by 
comparing the generated ones to the experimental data. 

Introduction 

Foraging Animals 

The characteristics of the movement of animals play a key role in a variety of ecologically relevant 
processes, from foraging and group behaviour of animals [1] to dispersal [2l[3] and territoriality [4]. 
Studying the behaviour of animals, simple random walk models have been proven effective in describing 
irregular paths [S]. While the first studies on random paths of organisms focused on uncorrelated step 
sequences |6], in many cases of studies of animal behaviour the directional persistence of the animals 
suggested a modelling in terms of correlated random walks (CRWs) OH]. In many complex environments 
an intermittent behaviour of animals is observed. In these cases an animal switches either randomly or in 
reaction to its environment between different movement patterns. The mechanisms which generate, and 
the factors which influence this switching behaviour have been shown to be important in understanding 
and modelling complicated animal paths [9l412j. While there is a source of switching between free flight 
and food inspections in the experiment we analyse |13| . here we concentrate on the former as detailed 
below. With no clear indication of additional intermittency, we will focus on non-intermittent models in 
the following. 



CRW/Reorientation Model 

The planar horizontal movement of an animal is often approximated by a sequence of steps: an angle 
a{t) describes the current direction of movement in a fixed coordinate frame, while the step length l{t) 
determines the distance travelled during a time step. The direction a{t), often determined by a specific 
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front direction of the animal, changes each time step by a random turning angle /3(t). The description of 
the dynamics in a co- moving frame, i.e., via the turning angle, turned out to be most useful for analysis of 
persistent animal movement [7]1H]. In many cases /3(t) is drawn independently and identically distributed 
(i.i.d.) from a wrapped normal distribution or a von Mises distribution |141ll5j for each time step, giving 
rise to a persistence in direction depending on how strongly the distribution is concentrated around 0. 
Usually the step length is taken to be either constant or it is drawn i.i.d. from some distribution. The 
step length can either be the result of a constant speed and a variable time step or (as in our case below) 
of a constant time step At and a variable speed s(t). 

This class of models can generate a variety of different dynamics. Two special cases with a uniform 
distribution for /? and a fixed time step At are the standard Gaussian random walk for step lengths 
l(t) = \z{t)\ where z is normally distributed and Levy flights for power-law tails in the step lengths 
distributions {l{t) ~ for 1 < < 3 and / > Iq). Related to Levy Flights, but using a time step 
proportional to the step length, are Levy Walks, which have been of interest as candidates for optimal 
search behaviour of foraging animals. They have been studied analytically [T^, by simulations [TU1I17) . 
and many experimental data sets have been statistically analysed to determine whether Levy Walks are 
suitable to describe the movement of animals (see, e.g., [I8H22] V 

As Levy-type models show anomalous diffusive behaviour, in contrast to models with a finite variance 
of the step length distribution and a fixed time step At, only the latter are included in the definition of 
correlated random walks which are also called reorientation models in the context of animal movement. 
Apart from pathological cases, CRWs are diffusive in the long time limit according to the central limit 
theorem. 

The estimation of the tortuosity of a trajectory is intimately connected to the distributions of the 
turning angle and speed [8lfT4l[23] . The relevance of the turning angle distribution for foraging efficiencies 
when searching in random environments has been analysed, e.g., in j24j . 

Generalisation of the Model 

In the following we will present a generalisation of the CRW above, which we then use to analyse 
bumblebee flight data. Given movement data with a constant time step At, the step length is determined 
by the speed s{t) = \v{t)\ of the animal. As we will be looking at a flying insect in a data recording using 
a small time step, we may expect to have a deterministic persistence due to the animals momentum. 
Additionally, the above CRW model assumes that s and /? are drawn i.i.d. which is sensible if At is large 
enough. However, for small time steps it cannot be excluded that the decision of the animal to turn left 
or right takes longer than the time step, which can correlate the turning angles /3(t) over a number of 
time steps. To allow for these possibilities we therefore model the changes in speed and turning angle via 
two coupled generalized Langevin equations. 



where we distinguish between the deterministic parts h and g and stochastic terms ip and (whose speed 
dependency will be discussed in the Results section). We assume that the noise processes are stationary 
with auto-correlation functions which may be non-trivial, and we make no further assumptions for the 
shape of their stationary distributions. While Eqs. (|l|2p represent a time-continuous description, the 
turning angle /3 still yields the change of a according to our fixed time resolution At. That is, /3(t) 
relates to a time-continuous angular velocity 7 of a via /3(t) ~ It-At l{'^)dT. The animals' position 
r(t) = (x(t), y(t)) is then given by dx/dt = scos(a(t)), dy/dt = ssin(a(t)) and da/dt = 7(t). Not having 
experimental access to 7, the numerical analysis is done with time-discrete data where the measured 
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Figure 1. Sketch of the foraging arena together with part of the flight trajectory of a 
single bumblebee. The bumblebees forage on a grid of artificial flowers on one wall of the box. While 
being on the landing platforms, the bumblebees have access to food supply. Of interest in this paper is 
the movement when the bumblebee is not near the flower wall. 

turning angle is given by f3{t) = a{t) - a{t-At) = Z{v{t),v{t - At)), where v{t) = {r{t + At) - r{t))/ At 
at times t = nAt, n £ N. 

Application to Experimental Data 

Analysing measured movement data of animals in their natural habitat is intricate due to a variety of 
factors which may influence the animal's behaviour, ranging from heterogeneous food source distributions 
[25H27] and predation threats [HIEH] to individual differences in behaviour within a population OH]. 
Here we analyse data obtained from a small scale laboratory experiment in which single bumblebees 
forage in an artiflcial flight arena |29j . The set-up is shown in Fig. [T] together with part of a typical 
trajectory of a bumblebee on its search for food. Each bumblebee can forage on an artificial flower carpet 
which is positioned on one of the walls of the arena. In this paper we are not interested in the behaviour 
resulting from the interaction with the flowers which has been studied in detail in |13j . Instead we only 
examine the search flights away from the flower carpet. (See section Materials and Methods for details.) 
We use our generalised stochastic model (Eqs. (|l|2p ) to describe these flights and to examine in which 
ways the behaviour deviates from a simple CRW model. Here we will focus on the horizontal movements. 
By neglecting the slower vertical movements, which are of more interest when analysing the starting and 
landing behaviour near flowers, we thus restrict ourselves to a two-dimensional model. 

Results and Discussion 
Estimation of Drift Terms 

Given the experimental data, we start determining the unknown parameters in our model by first estimat- 
ing the deterministic parts h{(3, s) and s) of the Langevin equation. This is done by numerical esti- 
mation [3OH33] of the components of the drift vector field (drift coefficients) D^{f3, s) = {g{l3, s), h{(3, s))^ 
of the corresponding Fokker-Planck equation via 
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Figure 2. Normalised drift vector field D^{f3, s) corresponding to the deterministic terms of the 
Langevin equations (Eqs. (|ll2p ) estimated via Eq. ([3]). The regular structure shows the quick relaxation 
to small angles, and the absence of strong cross-dependencies in the drift, i.e., the /3-dependence of the 
s-component of the vectors is weak and vice versa. 

where X = s)^ and (.) is the time average over the time series X conditioned on X{t) w X, where X 
is assumed to be stationary (for a detailed discussion see [32] )• The estimation of the drift terms is based 
on a Markov approximation: only those parts of the dynamics which match to a Markovian description 
in the state space variables f3 and s have their deterministic terms reflected in D^{X). Any other parts 
of the flight dynamics - stochastic as well as deterministic but not Markovian in (3 and s - are captured 
by the stochastic terms of Eqs. (|ll2p . Figure [2] shows the drift vector field, with normalised lengths of 
the vectors for better visibility. The nearly horizontal vectors show, that the drift quickly pushes the 
turning angle (3 towards 0, while the dynamics in the speed s is much slower. As the cross-dependencies 
of h{l3, s) on s and of <?(/?, s) on /? are weak, we can neglect them in our model. Since vector fields are 
hard to interpret, we will look at the projections in the following. 

Examining the drift h{/3) of the turning angle in Fig. |3] reveals that the drift term seems linear in j3 
— indeed we find numerically that its slope ~k matches exactly to a decay of the turning angle to in a 
single observation time step At by k w disregarding the noise term. This means that by integrating 
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Figure 3. Drift coefficient of turning angle. The deterministic drift /i(/3) as estimated from data 
(black) is in good approximation (red) linear in /?. (95% confidence intervals in grey) 
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Eq. ([T]) over a time At and approximating the drift h{f3) for small At by 
we have 



/3(t + At) - /3(t) = -k(3{t)/\t + 



t+At 



Ur)dT = -m + 



1+^' h{(i{T))dT ^ h{(i{t)W, 
t+At 



Ur)dT 



(4) 



With ^s(t) := Jt_^t £.s{T)dT and Eq. the time scale separation in the /3-Langevin equation due to the 
very fast relaxation means that we can simplify Eqs. (|ll2p to: 

$(t)=g(.(t))+^(t). 



(5) 
(6) 



While this reduction of the turning angle dynamics from dj3 / dt to (3 bears similarity to a simple reorien- 
tation model, the turning angles are still correlated and speed-dependent, as we will see below. 
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Figure 4. Drift coefficient of speed. The experimental deterministic drift coefficient g{s) (black, 
95% confidence intervals in grey) has been approximated by piecewise linear functions from one to three 
pieces (blue, green, cyan). The data shows the tendency to quickly increase low speeds. However, speeds 
above 0.27 m/s decrease more slowly, except for the rare high speeds. 

The speed drift g{s) displayed in Fig. |4] shows that the deterministic part of the speed-Langevin 
equation alone would have a stable fixed point around sq = 0.27 m/s. Comparing the slopes above and 
below So reveals that for s < sq the force towards sq is stronger than for s > sq. This is biologically 
plausible if one interprets sq as a preferred speed: if the bumblebee is slower it accelerates, but if it is 
faster it does not rush to decelerate as it would give up the energy spent to reach a high velocity. For 
very high velocities (over 0.55 m/s) the slope of g{s) increases again. This might be caused by the limited 
space available to the bumblebee in the flight arena. For our model we approximated g{s) by a piecewise 
linear function: 

5(s)«(s-so)x| ^^^^ , (7) 

where di > d2 > 0. As the very high velocities are rare, it made no difference in our model whether we 
used Eq. ([7]) or a piecewise linear function with three pieces. 

Velocity-Dependent Angle-Noise and Noise Auto-Correlations 

What we did not specify before was that the turning angle distribution may depend on the speed of the 
bumblebees. Given that the force a bumblebee can use to change directions is finite, the largest turning 
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Figure 5. Schematics of the dependence of /3 on speed s. Assuming a constant maximal force 
(circle) available to the bumblebee to accelerate during a time step, the distribution of the turning 
angle (3 depends on the previous speed St-i = \vt-i\- Illustrated is the change from large angles for low 
speeds (left) to a stronger concentration around 0° for higher speeds (right). 

angles have to be smaller when flying with high speeds (see Fig.[S|). This is consistent with the absence of 
simultaneously having high speed and large turning angle in the data - as is evident, e.g., from the data 
gaps in Fig. [2] However, animals can counteract this geometric dependence by varying the forces used 
for changing direction with the speed. We approximated the distribution for the turning angles for each 
speed s by a normal distribution. This approximation works best for low speeds. While there are some 
deviations for high speeds, it was not possible to reliably fit a better model due to the limited amount of 
data available. Figure [5] shows how its standard deviation depends on the current speed, ap decreases 
with increasing speed, however it does not decay to as a simple geometric model would predict (see 
Materials and Methods below). 




s [m/s] 

Figure 6. Log-log plot demonstrating the speed dependence of the turning angle 
distribution. The standard-deviation cr/j of the turning angle is shown as a function of the speed as 
estimated from data (black) and approximated by a shifted power-law (green) and a shifted exponential 
(blue). 95% confidence intervals for cr^ based on a x^-distribution are shown in grey. 

Instead cr^(s) decays roughly exponentially to a constant offset. We therefore model the turning 
angles as speed-dependent noise with a wrapped normal distribution [T3J[TS|: £,s{t) ^ Af{0,a^{s)) with 
cr^(s) = cie"'^^* -I-C3. This offset could cither be an effect of the boundedness of the flight arena, since the 
bumblebee has to turn more often to avoid walls when flying fast. Or it could be that the bumblebees use 
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stronger forces for turning during fast flights to maintain their manoeuvrabihty. It would be interesting 
to examine free-flight data to check for the cause. In other models in which the momentum of the animal 
is not important for the observed directional persistence, this cross-dependence is often neglected [7]. 




Delay time t [s] 

Figure 7. Log-log plot of auto-correlation of turning angles /?. The experimental data (black 
crosses) together with an exponential (magenta) and a power-law (blue) fit is shown with the large-lag 
standard error (grey). The green circles show the auto-correlation extracted from the simulated data. 

For the two stochastic parts of the Langevin equations, we estimated the normalised auto-correlation 
functions from the data. The turning angle auto-correlation is approximated by a steep power-law as 
seen in Fig. [71 which in this case is preferable to the alternative fit by a simple exponential decay. By 
subtraction of our approximation for the deterministic term g{s) from the observed speed changes ds/dt 
in Eq. ^ we estimated the distribution and auto-correlation of the noise term iplt) = ds{t)/dt — g{s{t)). 
In order not to overestimate the noise term, additive discretization errors of an approximate size of 
f error = Ax / At^ due to the finite resolution Ax = 10~^ m of the cameras have been accounted for, giving 
the variance cr^ = cr,^„oiay — cTorror- The noise term ip{t) is well approximated by Gaussian noise with 
an auto-correlation function acfl^'^(T) — ae^^^'^ + (1 ^ a)e^^^'^ (see Fig. [5]). While an auto-correlation 
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Figure 8. Auto-correlation of the non-deterministic speed changes i}){t). The auto-correlation 
function of ili{t) — ds{t)/dt — g{s{t)) estimated from the experimental data (dots) with two times the 
large- lag standard error (grey) and three fitted approximations: difference of 2 exponentials (red), 
difference of 2 power-laws (green) , difference of exponential and power-law (blue) . The outlier at 
T = 0.02 s is a discretization artifact due to the finite resolution of the data (see [34]). 
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function of the shape of acf^~^(T) = 6(r + 1)^p^ + (1 — 6)(t + 1)^^^ can be exluded, a difference between 
an exponential and a power-law acf^~^(r) ~ ce^'^^'^ + (1 — c)(r + 1)^^^ is not significantly worse than 
acfj^"*^. For our model we chose the simple difference of exponentials acf*^""^. 

As the observed anti-correlation between delays of 0.1s > t > 0.3 s happens on a time scale which is 
too short to be an effect of the boundedness of the experiment or of residual effects of the presence of the 
foraging wall |13] , it is unclear where the anti-correlation comes from. One could speculate that it might 
be the result of a stabilising mechanism in the bumblebee dynamics. 

Validation 

Given all the parameters of the full model (see Materials and Methods) estimated by minimizing the mean 
squared errors, we used them to generate artificial bumblebee trajectories, as follows: We simulated the 
dynamics using an Euler-Maruyama scheme with noise terms ^s(t),'0(i). In rare cases where the Gaussian 
noise ip{t) would lead to a negative speed despite the positive drift g{s) for s < sq, we enforce a non- 
negative speed by setting s{t) = 0. We correlated the noise terms in advance by modifying their power 
spectral density in the following way: we take uncorrelated noise of the wanted distribution, multiply its 
Fourier transform with the root of the desired power spectral density corresponding to our approximate 
auto-correlation function and then transform back |35j . To deal with the speed dependence of the turning 
angle noise S,s{t) we first correlate Gaussian noise and afterwards scale with (Ji3{s) at each time step in 
the integration scheme. While this does not reproduce the auto-correlation of the turning angle exactly, 
the error made is less than the errors from the estimation of acf^. A sample trajectory of a bumblebee 
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Figure 9. Simulated trajectory of a bumblebee. The complete model (Eqs. (|5l6p ) is simulated for 
200 s {— 10^ time steps) with an Euler-Maruyama scheme using already correlated noise for ^ and -0. 

simulated for 200 s using 10'"' time steps is shown in Fig. [HI Using the generated data we checked the 
validity of the model by comparison to the experimental data of all bumblebees. 

Figure ITOl compares the probability density function pdf(s) of the speed extracted from the simulated 
data with the experimental data. The auto-correlation functions of the speed and turning angle are shown 
in Figures [TT] and [T] Considering the number of rough approximations we have made for constructing 
our model, the agreement between simulation results and experimental data is very good. 

Summary 

We generalised a reorientation model which is often used to describe the correlated random walk of animals 
by explicitly modelling accelerations via Langevin equations. Analysing movement data from bumblebees, 
we extracted information on the deterministic and stochastic terms of Eqs. (|ll2p . Simulations of our model 
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Figure 10. Comparison of the speed distributions. The green (dashed) Une shows the 
probabihty density pdf (s) extracted from the simulated data, the black (solid) line shows the 
experimental data of all bumblebees (« 45000 data points). 
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Figure 11. Auto-correlation of bumblebee speed. The green (dashed) line shows the 
auto-correlation extracted from the simulated data, the black (solid) line from the experimental data 
with two times the large-lag standard error (grey). 

and comparison to the data have shown that the resulting model agrees very well with the experimental 
data despite the approximations we made for the model. With the estimation of the turning angle drift 
/i(/3) we found that while the usual assumption of i.i.d. turning angles is not valid in our case, the lack of 
a non-trivial drift and the weak auto-correlation of are consistent with the usual reorientation model. 
However, our generalised model exhibits significant differences in the non-trivial deterministic part g{s) 
of the speed change ds/dt and the speed dependence of the turning angles. In terms of active Brownian 
particle models |36[l37j we described the two-dimensional bumblebee movement by a particle with a non- 
linear friction term g(s) depending and acting only on the speed, driven by multiplicative coloured noise 
with different correlations for the angle component and the speed component of the velocity. While this 
combination of complications might make it difficult to treat the system analytically, progress into this 
direction has been made [351|35]. We remark that one could ignore the fast decaying auto-correlations 
of and 'ip{t) if one is not interested in the dynamics for short times, thus simplifying the model by 
using uncorrelated noise terms, since the effect of the noise autocorrelations on the long time dynamics 
is negligible. 

Given that the experiment which yielded our data is rather small and provided the bumblebees with an 
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artificial environment, it would be interesting to apply our new model to free-flying bumblebees to reveal 
how much the results depend on the specific set-up. This would clarify whether the fiight behaviour seen 
in the laboratory experiment survives as a flight mode for foraging in a patch of flowers in an intermittent 
model, with an additional flight mode for long flights between flower patches. The analysis of data from 
other flying insects and birds by using our model could be interesting in order to examine whether the 
piecewise linear nature of the speed drift and the trivial drift of the turning angle are a common feature. 
In view of understanding the small-scale bio- mechanical origin of flight dynamics, our model might serve 
as a reference point for any more detailed dynamical modelling. That is, we would expect that any more 
microscopic model should reproduce our dynamics after a suitable coarse graining over relevant degrees 
of freedom. 

Materials and Methods 
Experimental Data 

In this experiment 30 bumblebees (Bombus terrestris) were trained to forage individually in a roughly 
cubical flight arena with an approximate side length of 0.75 m. Figure [T] shows a diagram of the arena 
together with data from a typical flight path of a bumblebee. The flight arena included a 4 x 4 grid of 
artificial flowers on one of the walls. Each of the 16 flowers consisted of a landing platform, a yellow 
square floral marker and a replenishing food source where syrup was offered. For the analysis presented 
in this paper all data in zones (7 cm x 9 cm x 9 cm) around the flowers has been removed in order to 
analyse the search behaviour while foraging excluding the interaction with food sources. The 3D flight 
trajectories of the bumblebees were tracked by two cameras with a temporal resolution of At ~ 0.02 s. 
Each bumblebee was approximated as a point mass with a spatial resolution of 0.1 cm: its position was 
estimated by the arithmetic mean of all image pixels corresponding to the bumblebee via background 
subtraction. In total « 49000 data points were used for the analysis. For individual bumblebees an 
average of 51 search trajectories between flower zones have been sampled and analysed. The thorax 
widths of the bumblebees have a mean of 5.6 mm and a standard deviation of 0.4 mm. 

For calculating auto-correlations small gaps in the time series have been interpolated linearly. As the 
number of gaps was small the correlations for short times were not affected, however, the interpolation 
increased the usable data for long time delays. Trajectories were split at larger gaps, e.g., when entering 
a flower zone, to exclude correlations induced by flower visits. 

For a discussion of the influence of the boundedness of the flight arena and for the analysis of the 
foraging dynamics under varying environmental conditions see |13j . More details on the experimental 
setup can be found in [29 l l40 ) . 

Estimated Model Parameters 

The full set of parameters estimated from the data set which was used for the simulation is given here. 
For the deterministic drift of the speed the change of slope is at so = 0.275 m/s while the slopes are 
di ~ 0.16 and d2 = 0.06. The parameters for the standard deviation ct^(s) of the angle noise are 
ci = 126°, C2 = 12s/m, C3 = 12.5° and its auto-correlation is given by acf/3(T) — {t + i)-1-5476_ rpj^g non- 
deterministic changes 1^(1) of the speed are assumed to be normally distributed with standard deviation 
a^f, = 3.52 m/s^ and auto-correlated according to acf^~'^(r) where a = 1.44, Ai = 25.5 and A2 = 10.7. 

Speed Dependence of Turning Angles 

A simple model showing a dependence of the turning angles on the speed (see Fig. [5]) is given in the 
following. Assume that the velocity of an animal changes at each time step At by an acceleration vector 
which is given by a binormal i.i.d. random vector with variance cr^ in both directions. The turning angle 
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/3 between vt and Vt+At then depends on the quotient rjt := st/{V2a) between the former speed st — \vt\ 
and the noise strength a. By ehanging to the comoving frame of the animal and integrating out st+At 
the distribution of the turning angle is given by: 

p(/3) = + ^ — ?7C0s(/3)(l +erf(7?eos(^))) 

for — TT < /3 < TT. With vanishing speed s{t) = r](t) = the first term gives a uniform distribution as 
expected, and for ri{t) — > oo the distribution sharply peaks at (3 = with its variance cr^ approaching 0, 
similar to the behaviour of the simpler von Mises distribution. As the experimental bumblebee data does 
not show a decay to CT/j = but to a finite value (see Fig. [5]), this simple model does not hold: therefore 
the accelerations have to be modelled as speed-dependent. 
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